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ABSTRACT 


We present a comprehensive analysis of the Hubble Space Telescope observations of the atmosphere 
of WASP-121 b, a ultra-hot Jupiter. After reducing the transit, eclipse, and phase-curve observations 
with a uniform methodology and addressing the biases from instrument systematics, sophisticated 
atmospheric retrievals are used to extract robust constraints on the thermal structure, chemistry, and 
cloud properties of the atmosphere. Our analysis shows that the observations are consistent with a 
strong thermal inversion beginning at —10^ Pa on the dayside, solar to subsolar metallicity Z (i.e., 
—0.77 < log(Z) < 0.05), and super-solar C/O ratio (i.e., 0.59 < C/O < 0.87). More importantly, 
utilizing the high signal-to-noise ratio and repeated observations of the planet, we identify the following 
unambiguous time-varying signals in the data: 7) a shift of the putative hotspot offset between the 
two phase-curves and ii) varying spectral signatures in the transits and eclipses. By simulating the 
global dynamics of WASP-121 b atmosphere at high-resolution, we show that the identified signals are 
consistent with quasi-periodic weather patterns, hence atmospheric variability, with signatures at the 
level probed by the observations (~5% to ~10%) that change on a timescale of ~5 planet days; in 
the simulations, the weather patterns arise from the formation and movement of storms and fronts, 
causing hot (as well as cold) patches of atmosphere to deform, separate, and mix in time. 
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1. INTRODUCTION 


Spectroscopic measurements of transiting exoplanets 
have provided a wealth of “snapshot” information about 
the thermal structure, chemistry, and cloud properties 
of exoplanet atmospheres (e.g., Charbonneau et al. 2002; 
Tinetti et al. 2007; Kreidberg et al. 2014a,b; Swain 
et al. 2008b,a; Stevenson et al. 2014; Madhusudhan 
et al. 2014; Line et al. 2016; Sing et al. 2016; Tsiaras 
et al. 2018; Tsiaras et al. 2019; Benneke et al. 2019; 
Mansfield et al. 2021; Edwards et al. 2020; Skaf et al. 
2020; Pluriel et al. 2020; Mugnai et al. 2021; Line et al. 
2021; Changeat et al. 2022; Edwards et al. 2023; JWST 
Transiting Exoplanet Community Early Release Science 
Team et al. 2023). However, temporally varying infor- 
mation has yet to be unambiguously obtained by ob- 
servations. This is partly because, prior to the recently 
launched James Webb Space Telescope (JWST), exo- 
planet atmospheres have generally been studied with a 
single observation whose spectral feature signal-to-noise 
ratio (S/N) is too low. In an attempt to reduce the 
noise, the current standard practice is to average the 
signal from different observations; however, the averag- 
ing removes any temporal variability that may be cap- 
tured. On the other hand, when a single observation 
can achieve a high enough S/N, the observation of the 
planet is generally not repeated—due to observing time 
constraints. 

With the Spitzer telescope, a number of studies have 
analyzed repeated measurements of individual transiting 
exoplanets via photometric multi-epoch measurements 
of eclipses. Many of these studies did not detect at- 
mospheric variability below a certain level, due to the 
quality of the data (e.g., Agol et al. 2010; Crossfield et al. 
2012; Ingalls et al. 2016; Morello et al. 2016; Kilpatrick 
et al. 2020; Murphy et al. 2022). Others, however, have 
suggested time-dependent shifts in phase-curve offsets 
for at least three exoplanets—H AT-P-7 b, WASP-12 b, 
and Kepler-76 b—using either the Kepler or Spitzer tele- 
scopes (Armstrong et al. 2016; Bell & Cowan 2018; Jack- 
son et al. 2019; Wilson et al. 2021; Ouyang et al. 2023). 
The latter studies have speculated that such changes 
might be due to varying cloud structures, but conclu- 
sive interpretations of the datasets have remained elu- 
sive (Bell et al. 2019; Lally & Vanderburg 2022; Wong 
et al. 2022). Hence, presently there exists no unam- 
biguous detection of atmospheric variability in the at- 
mospheres of transiting exoplanets. 

In contrast, atmospheric variability is commonly re- 
ported for non-transiting exoplanets which are charac- 
terized by high-contrast imaging (e.g., Artigau et al. 
2009; Biller et al. 2015; Metchev et al. 2015; Biller 2017; 
Manjavacas et al. 2019; Vos et al. 2022). Among them, 


the ~11-19 Jupiter-mass planet VHS 1256-1257 b, which 
has recently been observed by the JWST-NirSpec and 
JWST-MIRI instruments as part of the Early Release 
Science program (Miles et al. 2023), exhibits one of the 
largest amplitudes of observed atmospheric variability; 
for example, a periodic brightness change of up to 3896 
with a period of —15 hours is reported by Zhou et al. 
(2022). Many other planetary mass companions exhibit 
similar levels of variability. 

Recently, intriguing possibility of time-variability for 
the ultra-hot Jupiter WASP-121 b has been reported in 
two studies, Wilson et al. (2021) and Ouyang et al. 
(2023). The former study compares spectra from a 
ground-based observation using Gemini-GMOS and a 
Hubble Space Telescope (HST) observation using HST- 
STIS and finds differences in the two spectra, which 
could be associated with a temporal variability. The 
latter study uses ground-based data from SOAR-GHTS 
and finds a spectra that also does not match that of 
the previous HST observations. These studies associate 
the observed differences with the presence of enhanced 
scattering slope in the case of GMOS, which could be 
explained by clouds or hazes, and varying abundances of 
molecular TiO/VO in the case of GHTS. Additionally, 
phase-curve observations with the HST (Mikal-Evans 
et al. 2022), Spitzer (Morello et al. 2023), and JWST- 
NirSpec G395H (Mikal-Evans et al. 2023) also report 
different phase-curve characteristics (i.e., “hotspot” off- 
set and shape), which could be indicative of moving 
hot regions in the atmosphere. However, the variabil- 
ity inferred in these works again relies on combining the 
constraints from different instruments and/or observing 
conditions. 

On the atmospheric dynamics modeling side, many 
hot Jupiter simulations in the past have suggested 
the presence of a single, stationary hot region east- 
ward of the substellar point—particularly between the 
~10*Pa to —10?Pa pressure levels (e.g., Cooper & 
Showman 2006; Dobbs-Dixon et al. 2010; Parmentier 
et al. 2018; Komacek & Showman 2020). However, at 
high-resolution, highly-dynamic, variably-shaped (and 
often multiple) hot regions emerge instead (Skinner & 
Cho 2021; Cho et al. 2021; Skinner & Cho 2022; Skin- 
ner et al. 2022). In these simulations, a long-lived giant 
storm-pair forms near the substellar point, drifts ini- 
tially toward one of the terminators, then rapidly trans- 
late westward thereafter—traversing the nightside and 
ultimately breaking up or dissipating near the eastern 
terminator; this cycle is quasi-periodic. "Throughout 
each cycle, hot (as well as cold) patches of air are chaoti- 
cally mixed over large areas and distances by the storms 
and sharp fronts around them. Similar mixing due to 
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storms and fronts has initially been predicted in high- 
resolution simulations (with a different initial condition 
than in the above studies) by Cho et al. (2003), who sug- 
gested that weather patterns would lead to potentially 
observable variability on hot exoplanets. 

In this backdrop, we present here results from an in- 
depth study of WASP-121b atmosphere—focusing on 
its variability. WASP-121 b is one of the best targets for 
atmospheric characterization because it is characterized 
by a high S/N and has been observed multiple times. 
It has been observed four times with the HST Wide 
Field Camera 3 Grism 141 (WFC3-G141): one transit 
in June 2016, one eclipse in November 2016, and two 
phase-curves in March 2018 and February 2019. Sig- 
nificantly, we utilize the entirety of these observations 
in this work. Previous analyses from a combination 
of facilities—HS T, TESS, Spitzer, JWST and ground- 
based facilities (Tsiaras et al. 2018; Evans et al. 2018; 
Mikal-Evans et al. 2020; Sing et al. 2019; Cabot et al. 
2020; Ben-Yami et al. 2020; Hoeijmakers et al. 2020; 
Borsa et al. 2021; Daylan et al. 2021; Mikal-Evans et al. 
2022; Changeat et al. 2022; Gibson et al. 2022; Azevedo 
Silva et al. 2022; Mikal-Evans et al. 2023)—have de- 
tected the presence of water vapour, absorbers of vis- 
ible radiation (VO and TiO), hydrogen ions (H^), and 
atomic species (Ba, Ca, Cr, Fe, H, K, Li, Mg, Na, V, 
and Sr); but, atmospheric variability has not yet been 
detected. 

The outline of this paper is as follows. Our basic 
methodology and codes are presented in Section 2. The 
reduction procedure for the construction of a consistent 
set of spectra from the raw observations is described in 
Section 3. Then, the results from our state-of-the-art 
atmospheric retrievals, which use the newly-developed 
“1.5D phase-curve retrieval" models (Changeat & Al- 
Refaie 2020; Changeat et al. 2021) and one-dimensional 
(1D) models, are presented in Section 4 and Section 5, 
respectively; the 1.5D models enable mapping of the at- 
mospheric properties (i.e. chemistry, cloud, and thermal 
structure) as a function of longitude using the entire 
phase-curve data, while the 1D models are used to ex- 
tract information from each of the recovered spectrum 
individually. In Section 6, we present the findings from 
the highest resolution, three-dimensional (3D) dynamics 
simulations of WASP-121 b atmosphere to date. Finally, 
discussion and conclusion are presented in Section 7. 
Additional material is included in an Appendix as well. 


2. METHODOLOGY 


'The WASP-121 b data we have analyzed in this work is 
obtained using HST WFC3-G141. This data is publicly 
available at the Mikulski Archive for Space Telescopes 


(MAST)!. Importantly, we have chosen to not include 
the observations from other telescopes or instruments, 
since the combination of multiple instruments is known 
to produce incompatible results (Yip et al. 2020, 2021; 
Edwards et al. 2023). The analyzed data is from one 
eclipse, one transit, and two phase-curves—comprising 
a total observing time of about 90 hours; each phase- 
curve observation contains two eclipses and one tran- 
sit. For consistency, we have used the same data re- 
duction pipeline, IRACLIS, and have adopted identical 
assumptions for all of the observations. The individual 
eclipse and transit events which are not from the phase- 
curve observations have been previously extracted us- 
ing IRACLIS by Changeat et al. (2022) [hereafter C22], 
and the phase-curve data has been recently analyzed by 
Mikal-Evans et al. (2022) [hereafter ME22]. However, 
the phase-curves have not been extracted using IRACLIS. 
'Therefore, we have re-analyzed the data with IRACLIS 
in order to ensure consistency of treatment with that by 
C22. 

Equipped with a consistently treated set of transit, 
eclipse, and phase-curve spectra, we use a suite of at- 
mospheric retrieval codes and a high-resolution atmo- 
spheric dynamics model code which is extensively tested 
and validated specifically for hot exoplanet simulations. 
Our aim here is to perform a robust extraction of the 
thermal structures and chemical abundance profiles in 
WASP-121 b's atmosphere, which allows us to estimate 
planetary formation markers and investigate potential 
time variability. Broadly, our methodology can be 
grouped into four main activities, or parts: 


Part 1: extracting a consistent set of WFC3-G141 light 
curves for the transit, eclipse, and phase-curve 
data using IRACLIS; fitting the light curves for 
the phase-curve data, using the POP (Pipeline 
of Pipes) code (see description in Materials and 
Methods); and, testing various assumptions to 
model the instrument systematics in the literature. 


Part 2: analyzing the recovered phase-curve dataset 
with the 1.5D retrievals developed by Changeat & 
Al-Refaie (2020) and Changeat et al. (2021) in the 
TAUREX3.1 framework (Al-Refaie et al. 2022a, 
2021), permitting time-independent, global prop- 
erties (e.g., mean metallicity Z and C/O ratio as 
well as mean thermal profiles at different longi- 
tudes) to be extracted. 


1 https://archive.stsci.edu/ ; 
data from HST proposals P14468 (1 transit, 02/06/2016, 
PI: Evans), P14767 (1 eclipse, 10/11/2016, PI: Sing), and 
P15134 (2 phase-curves, 12/03/2018 and 03/02/2019, 
PI: Evans) 


4 CHANGEAT, SKINNER ET AL. 2023 


—— 2018 visit |—— 2019 visit 
1.0015 
x »^ 
z $ lt E — 
v 1.0010 4 Be A 
D f : 
N e 
t 
5 1.0005 
a 
1.0000 


—0.6 —0.4 —0.2 
Time from mi-transit (days) 


0.0 0.6 


Figure 1. Corrected white light curves when the two WASP-121 b phase-curve visits are fitted individually. The two obser- 
vations, while close, present differences that could come from either atmospheric variability or instrument systematics. Note 
that the second-order coefficients from the sine phase-curve models are fixed to the best-fit values from the combined fit (see 


Figure B3). 


Part 3: analyzing the transit and eclipse data using 1D 
retrievals; and, incorporating the constraints (e.g., 
chemical parameters) obtained in Part 2 to reduce 
the degeneracies between temperature and chem- 
istry so that observations can be analyzed individ- 
ually rather than just in sum. 


Part 4: performing high-resolution, global atmospheric 
dynamics simulations with the psudospectral code 
BoB (Built On Beowolf) (e.g., Skinner & Cho 
2021; Polichtchouk et al. 2014), suitably optimized 
and forced with T-p profiles obtained in Part 2; 
and, interpreting the observed variability informed 
by these simulations. 


A more detailed description of each part is provided in 
the Materials and Methods section of the Appendix. 


3. SPECTRUM EXTRACTION FOR THE 
PHASE-CURVE DATA 


3.1. Combined white light curve correction 


Performing the extraction from the raw full-frame im- 
ages of the two observed phase-curves with IRACLIS, we 
obtain very similar results to ME22 (a comparison is 
provided in Figure A2). We then correct and fit the re- 
duced light curves with PoP (see description in Materi- 
als and Methods). Here, the two observations are fitted 
together, sharing orbital (mid-transit time tmia, inclina- 
tion i, and semi-major axis a) and model (planet-to-star 
radius ratios R,/R,, and phase-coefficients: Co, C1, C5, 
C3, C4) parameters, as well as the parameters for the 
short-term HST systematics. The parameters for the 
long-term HST systematics are not shared to accommo- 
date for the six different observation segments. Fitting 


the white light curves, we explore the effects of different 
short and long-term HST systematics on our results. 

For the short-term ramps, we have attempted two 
models: a simple exponential (e.g., as in Tsiaras et al. 
(2016b) [hereafter 2-param ISshort]) and a double ex- 
ponential (e.g., as in de Wit et al. (2018) and Mikal- 
Evans et al. (2022) [hereafter 4-param ISshort]). For the 
long-term systematics, three options are possible: lin- 
ear, quadratic, and hybrid (e.g., quadratic for the first 
segments of each visit and linear for the others, as in 
ME22). The comparison of these runs can be found in 
Figures B1 and B2. Overall, we conclude that assuming 
simple or double exponential short-term ramp does not 
change our results for this dataset. For consistency with 
ME22, we therefore adopt the double exponential ramp 
model for the remainder of the study; for the long-term 
ramp, however, the corrected phase-curve observations 
depend on the assumption of linear, hybrid, or quadratic 
option. 

Comparing the Bayesian evidence (E), a linear, hy- 
brid, and quadratic correction give log-values, In(E), of 
5733, 5797, and 5803, respectively. While the quadratic 
case gives a slightly higher log-evidence, such an assump- 
tion is too flexible for the second segment (i.e., the data 
around transit), generating artificial nightside flux and 
causing large degeneracies (see posterior distribution in 
Figure B2). We therefore adopt the more conservative 
approach and adopt a hybrid long-term ramp, as was 
done in ME22. The final recovered white light curve 
is corrected from the instrument systematics and com- 
pared with the results of ME22 in Figure B3. The resid- 
uals, in particular, demonstrate good agreement with 
the previous literature results. 
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3.2. Individual white light curve corrections 


To investigate the variability of the atmosphere and 
the instrument systematics in the available phase-curves 
observations, we have reproduced our white light curve 
analysis for each visit individually. Due to the lower 
amount of information, we choose to fix the orbital and 
second-order phase-curve parameters (C3 and C4) to the 
ones of the combined fit. For those runs, we have also 
employed the hybrid trend for the long-term ramps and 
the double exponential model for the detector short- 
term systematics. Figure 1 shows the corrected white 
light curves for the individual fits; Figure B4 shows the 
corresponding posterior distributions. 

'The recovered phase-curves in Figure 1 exhibit clear 
differences. For instance, the first-order phase of the si- 
nusoidal model is larger in the 2018 visit (0.19 rad) than 
in the 2019 visit (0.03 rad), or a variation in white tran- 
sit depth is found. These differences could come from a 
combination of atmospheric variability (e.g., movement 
of hot/cold regions or changing cloud coverage)—a pos- 
sibility we explore further in Section 6—and instrument 
systematic effects. In any case, these results suggest 
that combining HST phase-curve observations may not 
be straightforward. 


3.3. Spectral light curve fitting 


The spectral extraction is done using two different 
binning strategies, at “low” (i.e, as in ME22) and 
"medium" (i.e., as in C22) resolutions; see Materials 
and Methods for more details. Figures B5 and B6 show 
the corrected spectral light curves for the two cases, re- 
spectively. Inspecting the light curves, both strategies 
lead to similar residuals. We therefore moved on to the 
extraction of the spectra from the corrected light curves. 
Using 16 temporal bins (about 1.5h of observation), the 
final extracted spectra compared to that of ME22 are 
shown in Figure B7. Overall, the spectra agree very 
well—except for phase 0.07, where the reductions of our 
paper are consistent with each other but show larger 
flux than for phase 0.05 of ME22; this is despite the 
similar flux for phase 0.93, when compared with phase 
0.95 of ME22. This slight difference could arise from the 
fact that ME22 includes the in-transit planetary flux (af- 
ter removal of the transit signal) for bins 0.07 and 0.93 
(those bins labeled 0.05 and 0.95 in ME22 span 3h), 
while we chose to ignore the in-transit planetary flux for 
consistency and simplicity. 


4. PHASE-CURVE ATMOSPHERIC RETRIEVALS 


One of the goals of this study is to robustly charac- 
terize the bulk properties of WASP-121b atmosphere 
using the combined phase-curve data. As described in 


Materials and Methods, we use the 1.5D atmospheric 
retrievals to interpret the observed data.? This retrieval 
technique analyzes the two phase-curve observations us- 
ing a single unified atmospheric model (e.g., a single like- 
lihood); hence, it efficiently exploits all the information 
content—as opposed to the more traditional 1D retrieval 
performed individually for each phase (see, e.g., Steven- 
son et al. 2017). Since the “hotspot” offset Dys and (an- 
gular) size Agg are difficult to constrain from HST data 
alone (Changeat et al. 2021), we have tested different 
combinations and found that (Dys, Aus) = (30°, 50°) 
leads to the highest Bayesian evidence. Therefore, we 
here focus on this case. 

Figure 2 (and the associated animation) shows the 
best-fit spectra and recovered thermal structure from 
1.5D retrievals of the low-resolution data. Here, we 
obtain consistent 7'—p profiles by reproducing the re- 
trievals on the *Medium" resolution spectra as well as 
those from ME22 (see Figures C1 and C2), demonstrat- 
ing that this information is independent of the data re- 
duction. The chemical parameters (i.e., metallicity and 
C/O ratio) are slightly dependent on the spectral res- 
olution (especially in terms of precision, see Figure 3), 
which could be due to reduced correlation between the 
spectral channels at lower resolution, or information di- 
lution occurring during the fit due to the lower S/N of 
the light curves at higher resolution. Since our “Low” 
resolution reduction is consistent with ME22, we focus 
the rest of our discussion on this case; nevertheless, full 
posterior distributions for all three retrievals are pro- 
vided in Figure C3. 

Importantly, due to the resolving power of phase-curve 
data, our retrievals allow precise thermal and chemi- 
cal estimates at different locations in the planet’s at- 
mospheres to be obtained. Importantly, for ultra-hot 
Jupiters, the presence of absorption features for wa- 
ter, refractory species (TiO, VO, and FeH), and hy- 
drogen ions in WFC3 allows to break the degeneracies 
between metallicity and C/O ratio (see also Changeat 
2022). Given our retrieval assumptions, we find a strong 
thermal inversion on the dayside with the hottest region 
(here labeled hotspot in Figure 2) being ~300K hotter 
than the rest of the dayside between 10? Pa and 10? Pa 
(see red and orange profiles). 

'The thermal inversion could be caused by two differ- 
ent mechanisms. One mechanism is the production of 
energy at high altitudes by the presence of refractory 
molecules and H~ (the latter from Hz thermal dissocia- 


? See Changeat & Al-Refaie (2020), Changeat et al. (2021), 
and Changeat (2022) for additional examples of the 1.5D 
method. 
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Figure 2. Recovered temperature—pressure (T—p) profiles (left) and best-fit spectra (right) for the phases from 0.05 (blue) to 
0.5 (red), obtained from the phase-curve atmospheric retrieval. In the T—p plot, the shaded regions correspond to one and three 
sigma confidence regions (dark to light, respectively). The radiative contribution function is also shown in dashed line, colored 
for each region: hotspot (red), dayside (orange), and nightside (blue). These retrievals show good agreement with the observed 
data and demonstrate a strong dayside thermal inversion, with the presence of a hotter region (e.g. hotspot). The best-fit T — p 
profiles (solid lines, left) are used to thermally force the atmospheric dynamics simulations. This figure is accompanied by a 15s 
video, available online at the journal, showing the evolution of WASP-121 b emission (from blue to red) and the corresponding 
thermal structure as a function of phase. As the planet moves from transit to eclipse, absorption features in the data are 
replaced by emission features. These spectral variations enable the characterization of the thermal structure and chemistry 
across WASP-121 b's atmosphere. 
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Figure 3. Posterior distribution of the chemistry retrieved from the phase-curve data using the *Low" resolution spectra (Left), 
the *Medium" resolution spectra (Middle) and the spectra from ME22 (Right). The constraints obtained on metallicity Z and 
C/O ratio are consistent with a solar to slightly sub-solar metallicity and a super-solar C/O ratio, indicating the formation of 
the planet likely occurred beyond the snow-line. The blue line indicates the solar values for Z and C/O ratio. The full posterior 
distributions for the lou-resolution retrieval are available in Figure C3 
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tion): the temperature in the inversion region of the day- 
side is indeed hot enough to dissociate most molecules— 
including water and even more stable volatiles (CO and 
CO2) and refractory molecules (FeH, TiO, and VO), 
along with Hə (see the retrieved chemical profiles of 
Figure C4). At lower pressures, the atmosphere could 
partially be ionized, with an increased abundance of 
free electrons creating a continuum H^ opacity, as sug- 
gested for other similar ultra-hot Jupiters (Edwards 
et al. 2020; Pluriel et al. 2020; Changeat & Edwards 
2021; Changeat 2022). Another possible mechanism is 
heat deposition of breaking or saturating planetary and 
gravity waves launched from the atmospheric region be- 
low (e.g., Watkins & Cho 2010; Cho et al. 2015). Both 
mechanisms likely contribute to the observed thermal 
inversion layer. The retrievals we performed include a 
“gray” cloud model (i.e., constant opacity cloud deck); 
however, large cloud patches are not favored by the data 
(see Figure C3) despite the temperatures at the night- 
side being potentially suitable for silicate cloud forma- 
tion (Powell et al. 2018; Gao et al. 2021). 

Comparing the results from all the reductions (see also 
recovered T-p profiles and posteriors for other hotspot 
characteristics in Figures C5 and C6), we conclude that 
the data is consistent with a solar to slightly sub-solar 
metallicity Z and a super-solar C/O ratio. For the low- 
resolution case, we estimate the mean chemical charac- 
teristics of the planet to be log(Z) = —0.19*011 and 
C/O = 0.801993. A more conservative estimate, en- 
compassing the uncertainties from all the reductions and 
retrievals tested in this work is —0.77 « log(Z) « 0.05 
and 0.59 « C/O « 0.87. As a byproduct, this enables us 
to also speculate on the formation history of this planet; 
specifically, the obtained Z and the super-solar C/O ra- 
tio are suggestive of an early formation via significant 
gas accretion (i.e., without significant planetesimal pol- 
lution) and beyond the snowline of the proto-planetary 
disk (e.g., Óberg et al. 2011; Mordasini et al. 2016; Mad- 
husudhan et al. 2017; Brewer et al. 2017; Cridland et al. 
2019; Shibata et al. 2020; Turrini et al. 2021; Pacetti 
et al. 2022). 

'To further support the above conclusions, we have per- 
formed additional sensitivity tests, which are shown in 
Figure C7. We apply +30 departures, where ø is the 
retrieved uncertainty on the modified parameter, to the 
best-fit Z and C/O from the lou-resolution fit and com- 
pare the simulated spectra with the observations. In- 
troducing these departures leads to spectra that do not 
properly explain the observations, confirming the mag- 
nitude of the recovered uncertainties for the atmospheric 
chemistry of WASP-121 b. 


Due to the high constraints on the mean atmospheric 
properties of this planet obtained via the phase-curve 
data, the parameters extracted at this stage (e.g., ther- 
mal profiles and chemistry) can serve as priors for the 
subsequent parts of our analysis. In particular, assum- 
ing that Z and C/O ratio remains spatially homoge- 
neous and constant in time allows us to re-use the re- 
trieved values to analyze the transits and eclipse data 
individually; the thermal profile, chemistry, and cloud 
properties extracted from individual transit and eclipse 
observations are known to be much more degenerate 
on their own. Additionally, the recovered thermal pro- 
files provide important information for dynamics calcu- 
lations. For example, they can be re-introduced as an 
observation-driven forcing to enable physically realistic 
and case-specific simulations. 


5. TRANSIT AND ECLIPSE ATMOSPHERIC 
RETRIEVALS 


As mentioned, C22 has previously reduced the in- 
dividual transit and eclipse datasets with the IRACLIS 
pipeline; therefore, we make use of the spectra from 
that work directly. Already interesting differences ap- 
pear in the transit spectra, although the eclipse spec- 
tra look more alike (Figure 4). For both transits and 
eclipses, we perform 1D retrievals using the standard 
TAUREX3 models. We have first attempted to retrieve 
all the free parameters of the models without particular 
priors; but, as expected for HST data, the degeneracies 
between thermal structure, chemistry, and cloud prop- 
erties were difficult to break from individual HST tran- 
sit/eclipse spectra: we could not extract a consistent pic- 
ture. However, since Z and C/O ratio are expected to 
remain time-independent? and have better constraints 
from the phase-curve data, we have decided to re-inject 
this information from the 1.5D retrievals. Therefore, the 
chemistry is fixed to the median value from the retrieval 
on the low-resolution spectra; this has allowed to obtain 
a consistent fit of the spectra from all the visits (Fig- 
ure D1). For completeness, the posterior distributions 
are presented in Figures D2 and D3 for the transits and 
eclipses, respectively. 

For the transits, the three individual observations 
show the presence of covering hazes or clouds. Specifi- 
cally, the transit spectra captured during the two phase- 
curves (blue and green) are fully cloudy (i.e., consistent 
with featureless), while the observation obtained in 2016 


3 On considered timescales, the atmosphere is essentially a closed 
system. Note however that variable cloud formation via conden- 
sation can remove oxygen from the gas phase and locally change 
the C/O ratio. 
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Figure 4. Examples of observed transit (left) and eclipse (right) spectra of WASP-121 b analyzed in this work. The spectra 
shown are corrected for vertical offsets to highlight the differences in spectral shapes. Variability in the planet's weather patterns 
could create such variations in the spectroscopic data; this is strongly suggested in high-resolution simulations carried out for 
this planet in this work (see, e.g., Figures 6 and 7). For example, during the transits, the observed temporal variations could be 
interpreted as a formation of intermittent clouds and/or hazes; during the eclipse, the observed variations may be due to subtle 
changes in the thermal structure of the atmosphere—induced by motions of hot/cold regions from the planet's atmospheric 


dynamics. 


(orange) shows clear spectral modulations from water 
vapor. Note that the red end of this spectrum, however, 
cannot be fit properly using the chemical equilibrium as- 
sumption; however, free chemistry retrievals, which use 
H20, VO, and H^ opacities beyond equilibrium, could 
achieve a better match. 

Nevertheless, Transit 1 (orange) is not consistent with 
a featureless spectrum, as shown by A In(E) — 18.3, and 
possesses strong water vapor absorption features. Tran- 
sit 2 (blue) displays an interesting multi-modal solution. 
'The spectrum is best explained by either very high tem- 
peratures (T' zz 3000 K) at the terminator, forcing disso- 
ciation of the main molecules and leaving a flat contri- 
bution from the H^ continuum, or a more cloudy/hazy 
atmosphere with a slight slope towards the blue end 
of the spectrum. Given the un-physical nature of the 
high-temperature solution, we suggest that this second 
observation is consistent with the presence of hazes, es- 
pecially as a lower dimension featureless fit achieves a 
similar Bayesian evidence, Aln(E) — 0.5. For Tran- 
sit 3 (green), we find the observation consistent with a 
flat spectrum, which would be well explained by clouds 
and/or hazes, given Aln(E) = —0.6. In transit, stel- 
lar activity (i.e., unocculted stellar spots and faculae) 
can can cause important spectral variations between re- 
peated observations (Rackham et al. 2018; Thompson 
et al. 2023). However, long-term monitoring campaign 
from the ground (Delrez et al. 2016; Evans et al. 2018) 
suggests that WASP-121 is a very quiet and stable star. 
If confirmed, these results would indicate a transient 
formation of cloud/haze structures at the terminator of 
WASP-121 b. 


We note that gray clouds, which were only considered 
for the night-side, were not recovered by the retrievals 
on the phase-curve data. Many reasons could explain 
this result: i) the emitted signal in phase-curve does 
not probe the same altitudes as the transit data, ii) the 
observed clouds are poorly described by the gray cloud 
model, or iii) the clouds are located around the termi- 
nator region and not covering large patches of WASP- 
121b. Additionally, unambiguously determining the 
presence of gray clouds from emission data is more chal- 
lenging due to the degeneracies with the vertical tem- 
perature profile and the shorter geometrical path length 
through the atmosphere. 

For the eclipse data, the spectral differences are much 
smaller and difficult to infer by visual inspection of the 
spectra. We have conducted atmospheric retrievals and 
extracted the thermal structure for the five different 
eclipses individually (see left panel of Figure 5). The 
thermal profiles are overall consistent across the five 
observations (i.e., similar thermally inverted structure), 
but we find variations in the mean temperature of 310 K 
when averaging the profiles over the 10? Pa to 10? Pa 
region. Importantly, this range is much larger than the 
average one-sigma uncertainty of the profiles in the same 
region (the averaged standard deviation is 108K). For 
instance, in Figure 5, the thermal profiles extracted from 
Eclipse 2 (blue) and Eclipse 4 (red) are not consistent 
within the retrieved uncertainties. Similar to the phase- 
curve data, the observed differences in eclipse could be 
attributed to temporal variations of hot/cold regions in 
the planet’s dayside and/or a changing thermal struc- 
ture of the substellar point. 
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Figure 5. One-dimensional (1D) thermal structure recovered by our retrieval analysis of the five eclipse observations with one 
sigma confidence region (left), and T—p profiles from multiple times (t € [40,185] days) at the substellar point from a three- 
dimensional (3D) atmospheric dynamics simulation (right). The magnitude of variability in p € [10°, 10°] Pa is ~300 K, which 
is consistent with the variation predicted by the 3D simulation. Dashed gray lines show the vertical extent of the atmosphere 
modeled by the simulations in this study. Note, while these profiles are not like-for-like comparable because the retrieved thermal 
structure is global and substellar temperature predictions are local. 


While instrument systematics could remain (see Sec- 
tion above), the observed differences in hot region offset 
from the phase-curves, cloud coverage from the transits, 
as well as dayside thermal structure from the eclipses are 
all plausibly explained by the presence of atmospheric 
temporal variations. The observed differences are in fact 
expected from a theoretical atmospheric dynamics view- 
point due to the intense stellar heating contrast from the 
planet's parent star WASP-121. To investigate more 
precisely the possible origins of atmospheric temporal 
variations on the planet WASP-121 b and verify if they 
can affect our data to observable levels, we model its 
atmosphere with high-resolution dynamics simulations, 
to which we now turn. 


6. DYNAMICS MODELING 


We simulate the dynamics of WASP-121 b atmosphere 
with the BOB code at "T682L50" resolution, where 
T = 682 is the triangular truncation wavenumber (i.e., 
number of total and zonal modes each in the spherical 
harmonics) and L — 50 is the number of vertical layers 
(uniformly space in p); see Sec A.5, as well as Skinner & 
Cho (2021) and Polichtchouk et al. (2014), for detailed 
descriptions of the numerical model and simulation pa- 
rameters. The use of BOB at this resolution—to directly 
guide the retrieval interpretation with numerical robust- 


ness and verisimilitude—is a significant feature of this 
study. The simulations are performed to obtain a broad 
idea and insights into the variability plausible on planets 
like WASP-121b, when the flow is adequately resolved: 
it has recently been shown explicitly that hot-exoplanet 
simulations are not converged if the resolution employed 
is much below that in this work (Skinner & Cho 2021). 
As in most past studies, the atmosphere initially at rest 
is set in motion via a thermal relaxation to T—p profiles. 
Here the forcing profiles are obtained from the retrievals 
described above (e.g., Fig. 2) and prescribed. Profiles at 
many different times, which have deviated from the pre- 
scribed one due to the nonlinear atmospheric motion, 
are compiled in Fig. 5 (right panel) they should be 
compared with observations (left panel). 

Figure 6 shows temperature maps at p — 10? Pa from 
three widely separated times in the simulation. The 
maps demonstrate the highly variable nature of the 
planet’s temperature field at a pressure level from which 
observable flux would originate. Snapshots of the tem- 
perature at two different pressure levels, p = 5 x10? Pa 
and p — 10? Pa, show the vertical distribution, which is 
strongly barotropic—i.e, vertically aligned (Figure E1); 
a full movie of this simulation is provided in the Sup- 
plementary Materials. We also show chemical species 
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maps, which are simply post-processed using the instan- 
taneous density and temperature distributions from the 
simulation (Figures E2-E4) at p = 1 bar = 10? Pa and 
p = 50 mbar = 5x10? Pa (left two columns and right two 
columns, respectively), at t = 49 days and t = 62 days 
(left and right columns, respectively, for each p-level); 
distributions of the main relevant molecules (H5, H, 
H-, e^, H20, CO, CO2, CHy, TiO, VO, and FeH) are 
shown. Not surprisingly, variations (vertical, horizontal, 
and temporal) induced by WASP-121 b's atmospheric 
dynamics impacts the chemistry as well as temperature, 
the latter of which is discussed more in detail below. 
Note that here a simple chemical equilibrium assump- 
tion is made, which may not be valid everywhere in 
the modeled domain—particularly if the reaction time 
is comparable to the advection time. 

As can be seen in the figures, the complex motion of 
the atmosphere—and the organized structures therein— 
cause hot and cold regions to chaotically mix in time. 
Generally, the hottest region periodically forms slightly 
eastward of the substellar point, but it always moves 
away from the point of emergence. Depending on when 
the atmosphere is observed, the hottest region can even 
be located on the west side of the substellar point— 
either sequestered in a long-lived storm or generated 
near hyperbolic flow points between the storms (mid- 
dle and left panels in Figure E1, respectively). 

Interestingly, the dynamical behavior here is reminis- 
cent of that of WASP-96b in the p > 10* Pa vertical 
region. Skinner et al. (2022) have recently reported a 
quasi-periodic generation of giant cyclonic storms mov- 
ing away westward from the point of emergence. This is 
due to “deep heating” (i.e. strong heating at ~10° Pa), 
which may be experienced by some hot Jupiters. Here 
the similarity in behavior is likely due to the morpholog- 
ically similar T—p profiles in the aforementioned region 
on WASP-121b and WASP-96 b. The strength of the 
dayside-nightside difference is greater on WASP-121 b, 
but the profiles drive—and steer—the dynamics in a 
qualitatively similar way to WASP-96 b. 

'The similar dynamical behavior also leads to qual- 
itatively similar disk-averaged flux signatures at p — 
10? Pa, for example; cf. Figure 7 with Figure 4 in Skinner 
et al. (2022)*. In both figures, the fluxes exhibit quasi- 
periodic variations, with excess flux at the substellar 
(*Day side" in Figure 7) and east terminator regions. 
Hence, the flux is persistently “shifted to the east of the 
substellar point”—when averaged over the disk. It is 


^ NB., the former presents a black-body flux at 1.3m and the 
latter presents a much simpler flux (Sosa T^, where ogp is 
the Stefan-Boltzmann constant). 


important to note that, when not averaged, the hottest 
regions are rarely located near the equator and often sit- 
uated westward and at higher latitude of the substellar 
point (middle and left panels in Figure 6, respectively); 
the regions are not always vertically aligned either (cf. 
temperature maps from two pressure levels at Day 49 in 
Figure E2). Movies of the simulation show both clearly. 
'The movies also show that the timescale of the variabil- 
ity is ~5 planet days with sharp flux changes occurring in 
much shorter (< 0.5 day) windows—as was reported for 
WASP-96 b-like atmospheres by Skinner et al. (2022). 
'The above behavior overall may also explain why most 
infrared observations to date—except by Dang et al. 
(2018), Bell et al. (2019), and Morello et al. (2023)— 
have reported only *eastward-shifted hotspots". 

As expected, the magnitude of the variability in the 
model domain varies with the p-level: it is greatest at 
p ^ 10? Pa. Above p z 2x 10^ Pa, the dayside maintains 
a strong thermal inversion associated with an atmo- 
sphere that nevertheless remains highly variable. Given 
our model assumptions, this altitude-dependent behav- 
ior originates from the greater intensity of stellar irra- 
diation on WASP-121 b compared to that of a typical 
hot Jupiter, leading to a much shorter thermal relax- 
ation time, as well as the presence of the visible light 
absorbers H7, arising from the dissociation of H5 (Fig- 
ure E2). On WASP-121 b, while the stellar irradiation 
can still penetrate as deep as 10? Pa (as in a typical 
hot Jupiter), the upper layers of the atmosphere likely 
maintain much stronger dayside-nightside temperature 
gradient overall because of the shorter relaxation time 
(which is ox T~*: Andrews et al. 1987; Salby 1996; Cho 
et al. 2008). The thermal profile at the substellar point 
shows a short period (~5 days) temporal variation of the 
order of +200 K, which could be captured by the HST 
observations. 

To round out our investigation, we compare the tem- 
perature profiles and fluxes obtained in our dynam- 
ics simulations with the information obtained from the 
data. In Figure 5, we have already shown the agree- 
ment between the extracted temperature profiles from 
the observed eclipses (left panel) and the typical profiles 
resulting from simulations at the substellar point (right 
panel). Significantly, the predicted and observed spread 
of the temperature profiles at the substellar point in 
time is very similar, suggesting a qualitative agreement 
between observation and theory. From the simulations, 
we find that the three-o altitude-averaged temperature 
range is 680K between p € [10?,10?] Pa (e.g., assum- 
ing a normally distributed temporal spread of the T—p 
profiles, this value means that altitude averaged T-p 
profiles 340 K hotter or cooler than the mean should be 
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Figure 6. Instantaneous spatial temperature maps T'(A, $, p, t), where A is the longitude, ¢ is the latitude, p is the pressure, 
and t is the time); p = 10? Pa level at t — (49,62, 137) planetary days after the start of the simulation are shown. The maps 
demonstrate that very different temperature distributions arise at different times. Maps are in spherical orthographic projection, 
where the white circle marks the substellar point. The fields result from simulations which are thermally forced by the T—p 
profiles in Fig 2. The planet's *hotspot" is highly variable not only in location and time but also in shape. These storms are 
planetary-scale, coherent, and exhibit repetitive behaviors, leading to high amplitude and identifiable periodic signature in the 
planet's flux. 
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with variations of ~5 to ~10%. The normalized flux for the west terminator is generally lower than that of the east terminator, 
indicating that a phase-curve observation of this planet is more likely to have an eastward-shifted phase offsets. The level of 
variability is compatible with the uncertainties of our HST observations; however, note that the HST instrument systematics 
makes absolute flux measurements at such p-level highly uncertain (see, e.g., Yip et al. 2021; Edwards et al. 2023). 


considered outliers). In comparison, we find from the 
retrievals that the five eclipses have an average temper- 
ature range of 311 K, when averaged over the same pres- 
sure domain range. Despite the possibility of our data 
being affected by small HST systematics still remaining 
and the difficulty of directly comparing inherently differ- 
ent models (i.e, 3D vs 1D), the scale of the temperature 
variations from the eclipse observations and from the 
theory is compatible. 


As explained in the Materials and Method sec- 
tion of the Appendix, for three time frames, t — 
(49,62, 137} days, we have post-processed the simula- 
tion outputs to obtain the emitted spectral flux. Fig- 
ures E2—E4 demonstrate the wide range of physical and 
chemical conditions that could appear on WASP-121 b, 
which could strengthen the changes in cloud coverage 
and refractory chemistry claimed by Wilson et al. (2021) 
and Ouyang et al. (2023) respectively. In addition, we 


12 CHANGEAT, SKINNER ET AL. 2023 


present a comparison of the modeled spectra in Figure 8, 
which shows that variability of up to 1096 in the ob- 
served flux is compatible with the simulations. Such a 
level of flux variability agrees with theoretical works on 
other objects (Skinner et al. 2022), and is within the 
capabilities of HST—provided that its instrument sys- 
tematics are kept under control. Additionally, this level 
of variability should easily be captured by repeated ob- 
servations of similar planets from JWST and Ariel. 

Note that we obtain a lower level of variability when 
performing full radiative transfer calculation for the 
presented frames in Figure 8. This could be because 
of a selection bias (with the chosen frames) and/or be- 
cause of the baroclinicity (vertical non-alignment) of the 
atmosphere, which smooth the variability when the flux 
is integrated over a large pressure range. Moreover, we 
find that the variability is wavelength dependent, which 
we believe is caused by the changes in the chemistry. 
In particular, in the case of WASP-121 b, Ha thermally 
dissociates faster than H20; hence, we expect larger 
variability signals for wavelengths between 1.0 um and 
1.3 uum. 


7. DISCUSSION AND CONCLUSION 


The extreme atmospheric conditions of WASP-121b 
make it an ideal laboratory to test our understanding of 
physical and chemical processes in the atmospheres of 
exoplanets. Until now, detecting and studying weather 
patterns on exoplanets have remained elusive because 
of the lack of either adequate S/N or repeated, cross- 
verifiable observations that can be directly compared. 
However, multiple, comparable observations with the 
HST for WASP-121 b exist—and now permit progress 
to be made. Moreover, truly high-resolution, numeri- 
cally converged simulations also permit those observa- 
tions to be interpreted with some fidelity, since the gov- 
erning equations are now much more accurately solved 
than have been in the past. The uniform treatment 
in the analysis of observation, together with informed 
guidance from numerically accurate and validated sim- 
ulations, are the salient features of this study: without 
at least both, variability cannot be confidently assessed 
at present. Here, bolstered by state-of-the-art simula- 
tions, we identify spectroscopic variability in the HST 
observations of WASP-121 b. 

While some caution must still be exercised in inter- 
preting HST data, given the well-known high level of sys- 
tematics and the large number of assumptions required 
in reducing spectroscopic observations, we demonstrate 
here a strong potential evidence for variability associ- 
ated with weather on WASP-121b. Here the weather 
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Figure 8. Spectroscopic eclipse flux obtained at three differ- 
ent times, t = {49, 62, 137} days, when the dynamics calcula- 
tions are post-processed (top), and spectroscopic differences 
between two times, t = 62 days and t = 49 days. In the top 
panel, the eclipse spectra obtained for Observation 3 (green) 
and Observation 4 (purple) are also shown for reference. The 
simulation shows that the atmospheric variability should be 
wavelength dependent since, e.g., H^ is highly sensitive to 
thermal dissociation (impacting the short wavelengths) while 
H2O remains more chemically stable. 


is inferred from the following: i) the movement of the 
peak emission in two phase-curves, ii) the changing 
depth of the water feature in three transits, and i/i) the 
variable retrieved thermal profiles in five eclipses. On 
WASP-121b, the large dayside-nightside temperature 
gradient—which is not necessarily fixed in space and 
time?—is expected to power (as well as steer) dynam- 
ical, thermal, and chemical processes. These include 
vortex instability, gravity wave and front generation, 
thermal dissociation, chemistry changes, and potential 
cloud/haze formation (e.g., silicate clouds), whose con- 
sequences are observable. Other studies using ground- 
based data have also suggested variable atmospheric 
conditions on WASP-121 b (Wilson et al. 2021; Ouyang 
et al. 2023). 

A new finding in this study is that high-resolution dy- 
namics simulations forced by T(p) information retrieved 
from observations show that ultra-hot Jupiters, such as 
WASP-121 b, likely have hot regions that are generally 


5 due to the feedback from atmospheric dynamics or to the as yet 
incompletely understood thermal and/or orbital coupling with 


the host star 
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situated slightly eastward of the substellar point when 
disk-averaged—but whose actual shape and location are 
markedly variable in time. The variability is particularly 
visible in the modeled region, p € [10?, 10?] Pa, for which 
we have given the lion's share of focus in this study. 
This is in stark contrast with nearly all past hot Jupiter 
simulations, which show a fixed location and shape for 
a singular hot region; well-resolved simulations consis- 
tently indicate otherwise (e.g. Cho et al. 2003; Skinner 
& Cho 2021). More specifically, our simulations indicate 
that variability for WASP-121 b should be ~5% to ~10% 
in the disk-averaged flux with a frequency of ~5 days. 
Hence, the signatures generated by these quasi-periodic 
weather patterns should be detectable with current as 
well as future instruments—e.g., by the JWST (Greene 
et al. 2016) and Ariel (Tinetti et al. 2021) telescopes—if 
repeated, high-quality observations are obtained. 
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able from the Hubble Archive which is part of the Mikul- 
ski Archive for Space Telescopes. The specific observa- 
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A. MATERIALS AND METHODS 


A.l. Data reduction with IRACLIS 
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Figure A1. Displacement of the raw images onto the detec- 
tor for the 2019 observation. The three segments, separated 
by the re-acquisition events, are clearly visible on those di- 
agnostics. 


For the reduction and extraction of the spatially 
scanned spectroscopic images, we use the dedicated 
and publicly available pipeline IRACLIS (Tsiaras et al. 
2016b,c, 2018). The individual transit and eclipse events 
have already been extracted in the population study of 
C22, so we obtain those outputs from this study. For 
the phase-curves, however, the data has not previously 
been extracted with IRACLIS, so we perform all the ex- 
traction steps described in Tsiaras et al. (2019) as imple- 
mented in C22. Those steps consisted in zero-read sub- 
traction, reference-pixels correction, non-linearity cor- 
rection, dark current subtraction, gain conversion, sky 
background subtraction, calibration, flat-field correc- 
tion, and bad-pixels/cosmic-rays correction. We then 
use IRACLIS to extract the white (1.088-1.68 jum) and 
spectral light curves from the reduced images, taking 
into account the geometric distortions caused by the de- 
sign of the WFC3 tilted detector. The two observed 
phase-curves had re-acquisition events, separating the 
observations in three distinct segments (see Figure A1). 
Re-acquisitions cause displacements of the images onto 
the detector and induce larger systematics at the be- 
ginning of each event. Pixels with higher flux rates are 
affected more, causing wavelength-dependent long-term 
ramps that require an additional treatment. 

The extracted white light curves for both visits are 
compared to the ones obtained by ME22 in Figure A2. 
While completely independent, both reductions lead to 
very similar results. To extract the spectral light curves 


and study the impact of spectral binning on our conclu- 
sions, we explore two strategies: 


1) A medium spectral binning. 'To ensure consistency, 
we test the same binning as other IRACLIS reductions 
(Tsiaras et al. 2018; Changeat et al. 2022; Edwards 
et al. 2023), adopting an HST template composed of 18 
wavelength bins (see: Tsiaras et al. 2019). 


2) The low spectral binning. For better comparisons 
with ME22, we test their proposed 12 bins strategy. 


A.2. Light-curve analysis: the POP pipeline 


To fit the light curves, we have developed a new 
pipeline PoP! (Pipeline of Pipes), leveraging the flex- 
ibility of the TAUREX 3.1 framework (Al-Refaie et al. 
2022a). While TAUREX was originally designed for at- 
mospheric retrievals (Waldmann et al. 2015b,a), its last 
version provides generic classes that can easily be used to 
solve optimization problems outside of its original scope. 

In terms of code structure, PoP combines an obser- 
vation pipeline with a scientific model pipeline. Here, 
pipelines refer to a series of units representing individual 
transformation steps. In the context of this study, the 
"observation" pipeline has a unit to load the IRACLIS 
reduced light curves and a sequence of other units to cor- 
rect for the HST instrument systematics. The “model” 
pipeline contains an idealized light curve model. Note 
that the “model” pipeline can also possess additional 
transformation steps to convert the model’s outputs to 
the “observation” pipeline format. During optimization, 
the free parameters are marginalized over likelihood for 
both the observation and the model pipelines. This 
structure makes rapid modifications of the pipeline easy 
and creates a clear separation between astrophysical 
models and instrument models. 


Model Pipeline: In the case of WASP-121 b, we model 
the transit and eclipse events using the PYLIGHTCURVE 
package (Tsiaras et al. 2016a) as done in past works em- 
ploying IRACLIS, such as C22. The transit light curve 
drop (LC7) and the normalized eclipse light curve drop 
(LCg) from PYLIGHTCURVE are combined with a de- 
scription of phase-curve variations (LCs). Following 
the standard practice in the literature for HST (see e.g., 
Mikal-Evans et al. 2022), we model the phase-curve vari- 


10 Link to PoP: https://github.com/QuentChangeat/PoP. public 
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Figure A2. White light curves for the two WASP-121 b phase-curve visits obtained with our IRACLIS extraction and compared 
to the ones in ME22. The two visits were offset vertically (5x109 e^) and in time (by -257 orbital periods) for clarity. The 


observation re-acquisitions (Re-acq) are also annotated. 


ations using a first or second-order sinusoidal function: 


LCs = Cot c E cos(® cz) A [1 cos(2®—C4) |, 


(i) 
where ® is the orbital phase, Co is the minimum of the 
nightside flux, C1 is the maximum of the dayside flux, 
and C% is the phase-curve offset. C3 and C4 correspond 
to the second-order terms of the sine phase variations. 
As in Dang et al. (2018), the full phase-curve model 
Mec is constructed by: 


Mpc = LCr + LCs x LCg (2) 


In transit, we computed the limb-darkening coefficients 
with the Claret (2000) law using the stellar ATLAS 
models (Kurucz 1970; Howarth 2011; Espinoza & Jordan 
2015) included as part of the Exo TETHYS package 
(Morello et al. 2020a; Morello et al. 2020b). Those co- 
efficients were fixed during the light curve fits. 


Observation Pipeline: With HST, the instrument sys- 
tematics typically consists of a long-term trend affecting 
each visit and short-term ramps affecting each orbit. For 
the long-term trend (ISjong), the standard practice is to 
assume a linear or a quadratic behavior (Tsiaras et al. 
2016b; Kreidberg et al. 2018; Arcangeli et al. 2019). In 
the case of WASP-121 b, since re-acquisition events had 
occurred, we were required to fit the long-term trends 
separately for each segment of the observations (here 
segments are noted with the index ? with i € (1,2,3]). 
The corresponding long-term trend is: 


TSione = (1— Ab x (£t— 6) — Ad x (t — E) x N* (3 
g 0 s 1 s 


where Aj is the linear coefficient of the segment i, A1 the 
quadratic coefficient of the segment i, and N? is a nor- 
malization factor for the segment i. The time t, refers to 


the time at the beginning of each observation segment. 
For the detector short-term orbital ramps (ISghort), pre- 
vious studies have modeled its behavior using a combi- 
nation of exponentials. Most studies discard the first 
frame of each orbit due to the larger systematics. Addi- 
tionally, most methods also discard the first orbit, which 
is usually affected by a much larger and distinct expo- 
nential ramp. When discarding the first orbit and the 
first frame of each orbit, a simple exponential common 
to all the visits can be used (Kreidberg et al. 2014a; 
Stevenson et al. 2014; Tsiaras et al. 2016b): 


x) € 


1 


TSshort = 1+ Bo x exp ( 


where Bo and B, are the exponential factors, common 
to all visits in each observation and t? is the time at the 
beginning of each orbit ;. More recently, an evolution 
of this approach for the short-term ramps of HST has 
emerged (Zhou et al. 2017; de Wit et al. 2018). Their 
physically motivated model combines two exponential 
functions and is able to recover the first orbit. It is 


given by: 
t— t 
2 5 
a) 6 


TSshort = 1+ Bo x exp ( 


where R is the correcting function mainly acting on the 
first orbit and given by: 
t— to 
: 6 
Bs ) 6) 


By and B; are the exponential factors of R and t, is the 
time at the beginning of the visit. 


R= L+ Bo x exp ( 


In this paper, we explore the impact of various instru- 
ment systematics assumptions on the recovered data. 
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The approach, however, is the same in all cases and 
followed C22 for consistency. We start by fitting the 
white light curve with the goal to infer the free param- 
eters of the system that are wavelength independent — 
i.e., the mid-transit time (£54), the inclination (i), and 
the semi-major-axis (a). The other orbital parameters 
are fixed to literature values (Bourrier et al. 2019), as 
they have better constraints from complementary obser- 
vations. The fit is conducted with the nested sampling 
algorithm MULTINEST (Feroz et al. 2009; Buchner et al. 
2014), using 1024 live points and an evidence tolerance 
of 0.5. An initial estimate of the observational noise is 
made by IRACLIS (see Tsiaras et al. 2016b). However, 
we account for additional systematic noise or other un- 
accounted effects by re-scaling the uncertainties accord- 
ing to the RMS of the residuals, and perform a second 
fit using the updated uncertainties. In practice, this 
step has little effect for those WASP-121 b observations 
as the RMS of the residuals was very close to 1 at the 
first stage. Since nested sampling computes accurate 
Bayesian log-evidence, In(E), we note that the Bayes 
Factor can be used to perform model selection. 

To fit the spectral light curves, we employ a divide- 
white strategy similar to Kreidberg et al. (2014b); 
Tsiaras et al. (2016b); Mikal-Evans et al. (2022) and 
we divide each spectral light curve by the corresponding 
white light curve. Since the detector ramps are corre- 
lated between wavelengths, this step essentially removes 
the short-term ramps and reduces the baseline drift in 
the data. As such, the observation pipeline for the spec- 
tral light curve only contains the long-term trend cor- 
rection ISjong (ie., [Sshort is not needed). To match 
the divide-white observations, the model pipeline is also 
modified with an additional step. This step normal- 
izes the modeled spectral light curves, dividing by the 
median white light curve obtained during the prior fit. 
Each spectral light curve is then fitted separately with 
this model, keeping tmia, ? and a fixed to the best-fit 
value from the white. As with the white light curves, 
the errors are re-scaled to match the RMS of the resid- 
uals. 

Finally, we extract the phase-curve spectra from the 
binned corrected light curves using 16 different phase 
bins (4) of equal dimensions 0.05, giving us: ® € (0.07, 
0.12, 0.17, 0.23, 0.28, 0.32, 0.38, 0.43, 0.57, 0.62, 0.68, 
0.72, 0.78, 0.82, 0.88, 0.93}. Except for ® € {0.07, 0.93}, 
this matches the bins adopted in ME22. For those bins 
the planetary flux around ® = 0.0 is not included as it 
is blended during the transit, and a consistent bin size 
of 0.05 is used instead of 0.1 in ME22. Additionally, 
for this binning step, the finite integration time is not 
corrected for as the expected photometric distortions are 


below 5 ppm in the 0.05 intervals (Morello et al. 2022). 
For the transit spectra, self-blend (i.e, contamination 
by planetary emission) is not corrected for as the effect 
only affects the transit data by a few ppm in the case 
of WASP-121 b (Morello et al. 2019; Mikal-Evans et al. 
2022; Morello et al. 2021). 


A.3. Phase-curve retrievals 


We analyze the information contained in the phase- 
curve data using atmospheric retrievals. We employ the 
TAUREX 3.1 framework (Al-Refaie et al. 2021, 2022a), 
following the previously established methodology de- 
tailed in Changeat & Al-Refaie (2020); Changeat et al. 
(2021); Changeat (2022). We use the 1.5D phase-curve 
model, which is specifically designed to handle this type 
of observation, to simultaneously fit all the spectra in a 
Bayesian retrieval framework. The 1.5D model is com- 
posed of three different regions, referred to as hotspot, 
dayside, and nightside. Each region possesses indepen- 
dent properties allowing us to resolve large-scale atmo- 
spheric features from the data. The contribution of each 
region to the emitted flux at each phase is computed 
using a quadrature integration scheme (Changeat & Al- 
Refaie 2020). For this study, the structure of the planet 
is defined by 90 layers equally spaced in log space from 
p € [10°, 1071] Pa. To first-order, such a model accounts 
for the main one-dimensional biases discussed in Feng 
et al. (2016); Taylor et al. (2020); Changeat et al. (2021). 

As described in Changeat et al. (2021), the hotspot 
region is parameterized by two free parameters (hotspot 
size, Ays, and hotspot location, Dus). We have first 
attempted to retrieve those parameters but due to the 
large degeneracies between those parameters, the ther- 
mal structure, and the chemistry, this led to un-physical 
solutions (a similar behaviour was found in Changeat 
et al. 2021). Therefore, we have decided to fix those 
parameters and instead explored fixed values. For the 
hotspot size, we test Ays € {30,50} degree cases. For 
the hotspot location, we test Dys € (10, 20, 30, 40] de- 
gree eastward-shifted cases (those choices are also moti- 
vated by the results in ME22). 

To parameterize the thermal structure, the 
temperature-pressure (T—p) profiles are created by lin- 
early interpolating between T-p nodes (the pros and 
cons of such description are discussed in detail in the 
Appendix D of Changeat et al. 2021). For the hotspot 
region and dayside, we use seven nodes at fixed pres- 
sures (p € {10°,10°, 104, 103, 100, 10,0.1} Pa), while 
for the nightside, since the information content is re- 
duced due to the lower planetary emission, we choose 
to only use five nodes (p € {10°, 10°, 103, 10, 0.1} Pa). 
For the chemistry, we employ the GGCHEM (Woitke 
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et al. 2018) chemical equilibrium code in its TAUREX 
3.1 plugin and we couple the two only free parameters 
(metallicity and C/O ratio) between the three different 
regions. Previous works have shown the importance of 
dis-equilibrium processes for hot Jupiters (Moses et al. 
2011; Drummond et al. 2020; Venot et al. 2020, 2012; 
Al-Refaie et al. 2022b), however, given the temperatures 
of WASP-121 b, chemical reactions should be fast, fa- 
voring chemical equilibrium for the investigated species 
(Parmentier et al. 2018; Kitzmann et al. 2018). The 
radiative transfer includes absorption from the main 
expected opacity sources via ExoMol line lists (Ten- 
nyson & Yurchenko 2012; Chubb et al. 2021), namely: 
H20 (Polyansky et al. 2018), CO (Li et al. 2015), CO2 
(Yurchenko et al. 2020), CH4 (Yurchenko et al. 2017), 
TiO (McKemmish et al. 2019), VO (McKemmish et al. 
2016), FeH (Bernath 2020), and H^ (John 1988; Ed- 
wards et al. 2020). We also consider Collision Induced 
Absorption (CIA) by Ho-H» and H2-He pairs (Abel et al. 
2011, 2012; Fletcher et al. 2018) and Rayleigh scattering 
(Cox 2015). A fully opaque cloud deck (referred here as 
gray clouds) was also included on the night side of the 
planet but we were not able to find evidence for clouds 
from this phase-curve data. The parameter space of the 
model is explored using the nested sampling optimizer 
MULTINEST (Feroz et al. 2009; Buchner et al. 2014), 
with 500 live points and an evidence tolerance of 0.5. 
'To explore the parameter space, priors were chosen to 
be non-informative (i.e., uniform priors). Specifically, 
for the temperature points, we explored the space from 
T € [300, 6000] K. For the chemistry, the metallicity was 
explored in log space from Z € [0.1,100] times solar, 
while the C/O ratio was explored from C/O € [0.1, 2]. 
From this retrieval analysis, we were able to extract 
averaged chemical properties of WASP-121 b as well as 
the thermal structure of the three considered regions. 


A.4. Transit and eclipse retrievals 


To evaluate the variability of WASP-121b’s atmo- 
sphere, we also analyze each eclipse and transit spectra 
individually, using 1D atmospheric retrievals with TAU- 
REx 3.1. Previous studies have shown the difficulty 
of extracting reliable constraints from single HST vis- 
its (Changeat et al. 2020) due to degeneracies between 
chemical abundances and thermal properties. To break 
those degeneracies, we use our most accurate estimate 
of the time-independent parameters from our phase- 
curve retrieval as priors for our individual retrievals. 
Using the low-resolution results, the metallicity Z of 
the atmosphere is therefore fixed at log(Z) = -0.19, 
while the carbon-to-oxygen ratio C/O is fixed at C/O 
= 0.80. Note that this simplification is not expected to 


always be correct, for instance cloud condensation can 
locally (and temporally) change the C/O ratio of the 
gas phase, which we do not model here (i.e., GGCHEM 
is used without condensation). However, without ad- 
ditional knowledge or constraints on condensates in 
WASP-121 b, this remains a reasonable and necessary 
assumption. 


Transits: Since HST transits probe a narrow pres- 
sure range, and because it is mainly affected by the 
atmospheric scale height (Rocchetto et al. 2016), we 
consider a simple isothermal profile with a unique free 
parameter T for those observations (spectra shown in 
Figure D1). As with the phase-curve data, the tempera- 
ture is explored with uniform priors (T' € [300, 6000] K). 
However, transit is much more sensitive to clouds than 
eclipse or phase-curve, so we have used a more complex 
representation of clouds from Lee et al. (2013), for which 
particle size and mixing ratio were fitted. This cloud 
model was favored by the Bayesian evidence compared 
to the gray-cloud case. 


Eclipses: As planetary radius is known to be de- 
generate with temperatures in HST eclipses (Edwards 
et al. 2020; Pluriel et al. 2020), we fix this parameter 
to the literature value. For each spectrum (see Fig- 
ure D1), we retrieve a thermal profile with a similar 
parameterization to the phase-curve case. Namely, the 
T-p profile is parameterized by linearly interpolating 
between seven freely moving T-p nodes. The pressure 
of each node is fixed to log-spaced values of pressures 
(ie., p € (106,105, 104, 103, 100, 10, 0.1} Pa) and we re- 
trieve the temperature of each point individually using 
the same uniform non-informative priors. We refer to 
Changeat et al. (2021); Rowland et al. (2023) for a more 
complete discussion on thermal structure parameteriza- 
tions and their trade-offs. A simpler 3-point thermal 
profile (where the pressure levels of each node are left 
free) was also tested, which did not change our overall 
conclusions. For the eclipse spectra, we also decided to 
run a simple blackbody planet fit, which served as our 
comparison baseline. Following this procedure, and be- 
cause of the additional chemistry priors from our phase- 
curve analysis, we have obtained well-defined thermal 
structures for each of the five eclipses. 


A.5. Dynamics modeling 


We model the atmospheric dynamics of WASP-121 b 
using the pseudospectral dynamical core, BOB (e.g., 
Rivier et al. 2002; Scott et al. 2004; Polichtchouk et al. 
2014; Skinner & Cho 2021; Skinner & Cho 2022). The 
core has been outfitted and set up in Skinner & Cho 
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(2022), Skinner & Cho (2021), and Cho et al. (2021) es- 
pecially for high-resolution hot Jupiter simulations; and, 
we refer the reader to those works for a more complete 
description of the code and governing equations solved. 
However, for the readers’ convenience, we provide a brief 
summary of the key features of the model here — espe- 
cially as they pertain to modeling of WASP-121 b at- 
mosphere. Other numerical models have been used to 
study hot Jupiter atmospheric dynamics in the past (e.g. 
Showman & Guillot 2002; Cho et al. 2003; Cho et al. 
2008; Cho 2008; Dobbs-Dixon et al. 2010; Thrastarson 
& Cho 2010; Polichtchouk et al. 2014; Heng et al. 2014; 
Mayne et al. 2014; Mendonca et al. 2018; Parmentier 
et al. 2018, and references therein), and we direct the 
reader to consult those works for instructive context. 

BoB calculates the large-scale dynamics of the atmo- 
sphere by numerically solving the traditional and hydro- 
static primitive equations in (longitude, latitude, pres- 
sure) = (A, ¢,p) coordinates an in vorticity-divergence— 
potential-temperature formulation. In the vertical (p) 
direction, BoB employs a second-order finite differenc- 
ing scheme with free-slip boundary conditions at the 
top and bottom pressure surfaces. The present study 
builds upon extensive testing and validation of BoB, 
including simulations at high-resolution and under nu- 
merically challenging conditions resembling those found 
on WASP-121b (Polichtchouk & Cho 2012; Polichtchouk 
et al. 2014; Cho et al. 2015; Skinner & Cho 2021). 

For the physical setup of WASP-121b we use the 
parameters in Delrez et al. (2016). To simulate irra- 
diation from the planet's host star we implement an 
“idealized” thermal forcing using the Newtonian relax- 
ation scheme which accelerates the initially resting at- 
mosphere (u = 0) towards specified hotspot and night- 
side equilibrium 7—p profiles that are obtained from the 
phase-curve retrievals of the planet (see Fig. 2). The 
initial temperature is the average of the hot spot and 
nightside equilibrium temperatures. Due to the planet’s 
close proximity to its star, the effect of its spherical ge- 
ometry on stellar irradiation deposition is accounted for 
using a cosine profile to graduate the hot spot tempera- 
ture from the substellar point to the terminators. The 
radiative cooling time 7,(p) is computed from the initial 
temperature profile following Cho et al. (2008); 7;(p) is 
approximately linear in log(p), ranging from O(109)s at 
p = 10? Pa to O(107)s at p = 10? Pa. 

For the numerical setup, we use high horizontal and 
high vertical resolutions: T682 and L50 respectively. 
In the former, ‘T’ denotes the highest wavenumber re- 
tained in the spherical harmonic expansion (the trunca- 
tion wavenumber) is n, = 682 and the latter ‘L’ denotes 
the number of vertical levels which are distributed lin- 


early in p-coordinate. By “high-resolution” we mean 
our simulations are above the minimum resolution re- 
quired for numerical convergence of flow solutions on 
hot, tidally synchronized exoplanets (Skinner & Cho 
2021). Note that high vertical resolution is also nec- 
essary to ensure the baroclinic region of the dayside and 
hotspot temperature profiles (p < 10^ Pa) is well rep- 
resented in p. A small timestep size of At = 6 seconds 
is used concomitantly with the fine grid spacing to en- 
sure flows at the maximum sound speed (i.e., c, ~ 4880 
m/s in the hottest region of the atmosphere) are well 
captured. Hence, the simulations maintain a Courant- 
Friedrichs-Lewy (CFL: Courant et al. 1928) condition 
of well below unity. Simulations are time-integrated for 
200 planet days, to ensure they reach a quasi-stable state 
long after their initial ramp-up period of ~40 WASP- 
121 b days. 

For the domain size, we model the expected p range 
(from p, = 10? Pa to p = 10° Pa) from which radiation 
originates, as indicated by HST data (see Figure 2 and 
Figure 5). We have verified that the model equations 
are valid for this region by confirming that the retrieved 
thermal forcing profiles are stably stratified (i.e. that the 
Brunt-Váisálà frequency M = y (g/p) (dp/dz) is real). 
This is shown in Fig. C2. Below this (p > 10°), the 
nightside temperature profile exhibits a jump in strati- 
fication. Although this could be due to increased uncer- 
tainty on the retrieved thermal structure in this region. 

Finally, for numerical dissipation we use a high-order 
hyper-viscosity V!9 with a small corresponding artifi- 
cial viscosity coefficient of vig = 10478. This prevents 
excessive kinetic energy removal from small-scale flows 
and hence ensures the dynamics of large-scale flows are 
well represented (see e.g. Cho & Polvani 1996; Skinner 
& Cho 2021). In addition, a very weak Robert—Asselin 
filter with coefficient e = 0.02 is used to filter the ad- 
ditional computational mode arising from the models 
leapfrog time integration scheme (Asselin 1972; Thras- 
tarson & Cho 2010). Besides this, no other drags are 
applied as these can coerce the flows to dynamically un- 
physical states (Polichtchouk et al. 2014). The simula- 
tions are allowed to evolve freely under thermal forcing 
from the retrieved temperature profiles. 

For the post-processing of the three-dimensional at- 
mospheric dynamics simulations, we employ two differ- 
ent approaches for evaluating the evolution of WASP- 
121,b from an observational perspective. The first is a 
one-dimensional time-series analysis of flux emitted by a 
single layer at the mid-region of our computational do- 
main (0.1 bar) in order to evaluate the qualitative behav- 
ior of the planet's weather over hundreds of planet days. 
Here the black-body emission is calculated from BoB 
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temperature maps, centered on key regions of interest 
(i.e., the substellar point, eastern terminator, western 
terminator, and antistellar point). While this approach 
does not account for the entire domains contribution 
towards atmospheric variability, it enables the key dy- 
namical processes in the atmosphere to be isolated and 
studied in detail. 

For second phase of our post processing, we link the 
outputs of BOB with the TAUREX3 library to simulate 
observables and produce detailed 3D chemical maps of 
the planet's atmosphere. Due to the large size of the 
BoB simulation grid cube we isolate frames from the 
BoB calculations with significant spatial temperature 
differences that are likely observable for this analysis 
(e.g., t = 50 and t = 62 days). First, we derive the 
chemical maps over the entire (A, 9, p) = (2048, 1024, 50) 
(i.ee., grid cube) by performing calculations with the 
GGCHEM (Woitke et al. 2018) chemical equilibrium 


code using the T-p values from each grid square of 
the BoB simulations. We maintain fixed values for the 
metallicity and C/O ratio based on the median values 
obtained during our atmospheric retrieval exploration 
(log(Z) = -0.19 and C/O = 0.80). The flux emitted 
by the planet is then computed using the TAUREX3 
plane-parallel radiative transfer model, modified for our 
three-dimensional grid and to account for the changing 
viewing angles. That is, for each column of the compu- 
tational grid, the flux is propagated upwards from 0.5m 
to 504m at resolution R = 15,000 and then summed by 
taking into account the viewing angle of each grid ele- 
ment. In this calculation, the same opacities as during 
the retrievals were used: molecular absorption via EXO- 
Mor cross-sections, H^ opacity, Collision Induced Ab- 
sorption (CIA) and Rayleigh scattering. We computed 
the planetary flux in those frames as if the planet were 
observed in eclipse (i.e., phase 0.5). 
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B. COMPLEMENTARY FIGURES TO THE SECTION 3 


'This appendix contains the complementary figures to the main article Section 3, Figure B1 to Figure B7. 
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Figure B1. Comparison of recovered white light curves, fitted with the POP when using different assumptions for the detector 
and long-term ramp models. We show in solid line the median of each tested model, and for our preferred set of assumption, 
we show the resulting observations. Red: Two parameters exponential ramp model from Tsiaras et al. (2016b) and hybrid 
long-term ramp; Blue: Four parameters exponential ramp model from de Wit et al. (2018) and hybrid long-term ramp; Green: 
Two parameters exponential ramp model and linear long-term ramp for each segment; Orange: Two parameters exponential 
ramp model and quadratic long-term ramp for each segment. The hybrid long-term ramp consists in a quadratic trend for the 
first segment of each visit and a linear ramp for the remaining four segments, as done in ME22. This, as well as Figure B2, 
demonstrate that the recovered light curve depends on the assumption for the long-term ramp. 
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Figure B2. Posterior distributions of the white light curve fits with PoP using four different detector ramp models (see 
Figure B1). Red: Two parameter exponential ramp model from Tsiaras et al. (2016b); Blue: Four parameters exponential 
ramp model from de Wit et al. (2018). Green: Two parameters exponential ramp model and linear long-term ramp for each 
segment; Orange: Two parameters exponential ramp model and quadratic long-term ramp for each segment. The recovered 
orbital parameters are independent of the detector short-term ramp, but they are impacted by the choice of the long-term trend. 


ATMOSPHERIC VARIABILITY OF WASP-121B 27 


1.0015 4 $ 


1 Evans et al. 2022 
$4 2018 visit 


1.0010 4 ne 
4 2019 visit 


Normalized Flux 


Normalized Flux 


Residuals 


—0.6 -0.4 —0.2 0.0 0.2 0.4 0.6 
Time from mi-transit (days) 


Figure B3. Corrected white light curves when the two WASP-121 b phase-curve visits are fitted together. We show our results 
with the IRACLIS extraction and compared to the ones in Mikal-Evans et al. (2022). Top: Corrected light curves zoomed in the 
eclipses; Middle: Full corrected light curves; Bottom: Residuals between our best-fit light curve model and corrected data with 
red and blue for the IRACLIS reductions and green for the reduced data in Mikal-Evans et al. (2022). 
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Figure B4. Posterior distributions of the white light curve fits with POP when fitting the two phase-curve visits independently. 
''he recovered parameters show differences in both the phase-curve model and the instrument systematics. 
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Figure B5. Corrected spectral light curves (left) and corresponding residuals (right) for the low spectral binning. This is the 
same binning employed in Mikal-Evans et al. (2022). Higher binning resolution fits (e.g. medium resolution) are also performed 


and presented in Figure B6. 
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Figure B6. Corrected spectral light curves (left) and corresponding residuals (right) for the medium spectral binning. The 
binning is similar to previous works using reductions from IRACLIS, such as Tsiaras et al. (2018); Changeat et al. (2022); Edwards 


et al. (2023). 
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Figure B7. Emission spectra of WASP-121 b, obtained at different phases from various reduction methods. Green: spectra 
at low resolution from ME22. Red: spectra at low resolution obtained using the 4-short instrument systematic model. Blue: 
spectra at medium resolution obtained using the 4-short instrument systematic model. Orange: spectra at medium resolution 
obtained using the 2-short instrument systematic model. All the reductions are consistent with each other. Note that phases 
0.07 and 0.93 do not exist in ME22; hence, we instead plot their phase 0.05 and 0.095. 
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C. COMPLEMENTARY FIGURES TO THE SECTION 4 


This appendix contains the complementary figures to the main article Section 4, Figure C1 to Figure C7. 
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Figure C1. Temperature — pressure profiles (T' — p) obtained by the 1.5D retrievals on the spectra from different reductions. 
Left panel: nightside. Middle panel: dayside. Right panel: hotspot. We show the extent of the radiative contribution function 


for the low-resolution retrieval with dashed lines. The retrievals on those different reductions are consistent and provide a similar 
picture. 
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Figure C2. Profiles of Brunt-Váisálà frequency squared A? = (g/p) (dp/dz) for the retrieved thermal profiles in Fig.C1. 
Dotted lines at p = 10! and p = 10° Pa show the upper and lower boundary of the GCM model. In this region, flows are 
stably stratified (V? > 0); hence, they satisfy the hydrostatic balance approximation of the model equations. For p > 10°, the 
nightside, hotspot and initial profiles exhibit a jump in stratification. 
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Figure C3. Posterior distributions for the atmosphere of WASP-121 b obtained by the 1.5D retrievals on different reduction 
methods. Green: spectra at low resolution from ME22. Red: spectra at low resolution obtained using the 4-short instrument 
systematic model. Blue: spectra at medium resolution obtained using the 4-short instrument systematic model. 
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Figure C4. Chemistry of the main species recovered from the phase-curve data (low-resolution reduction). Red: hotspot; 
Orange: dayside; Blue: nightside. Those chemical profiles show thermal dissociation of main molecules such as H2, H20, CO, 


CO», TiO and VO, at the hotspot of WASP-121 b. 
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Figure C5. Temperature — pressure profiles (T — p) obtained by the 1.5D retrievals on the low-resolution spectra by varying 
the hotspot size (Ans) and offset (Dus). Left panel: nightside. Middle panel: dayside. Right panel: hotspot. While the 
thermal structure are different, the conclusions on the thermal structure of WASP-121 b from those runs would be the same, 
independently from the hotspot assumptions. We also note that one model (Ans = 50 and Dys = 30) has a significantly higher 
Bayesian evidence. The models obtained: In(E) = 1550.3 for Aus = 30 and Duys = 30, In(E) = 1510.6 for Aus = 50 and Dus 
= 10, In(E) = 1538.8 for Aus = 50 and Dus = 20, In(E) = 1554.5 for Aus = 50 and Dys = 30, and In(E) = 1537.5 for Aus = 
50 and Dus = 40. 
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Figure C6. Posterior distributions for the atmosphere of WASP-121b obtained on the low-resolution spectra by varying the 
hotspot size (Ans) and offset (Dus). The color code is the same as in Figure C5. Independently from the hotspot assumptions, 
all those retrievals of WASP-121b phase-curve data have a super-solar C/O with 0.62 < C/O < 1.11, while the retrieved 
metallicity is between C/O with -1.27 < log(Z) < 0.77. For all cases, we do not find evidence of a fully opaque nightside cloud 
deck. One model (Ans = 50 and Dys = 30) has a significantly higher Bayesian evidence. The models obtained: In(E) = 1550.3 
for Aus = 30 and Dys = 30, In(E) = 1510.6 for Aus = 50 and Dys = 10, In(E) = 1538.8 for Aus = 50 and Dug = 20, In(E) 


1554.5 for Aus = 50 and Dys = 30, and In(E) = 1537.5 for Aus = 50 and Dus = 40. 
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Figure C7. Sensitivity tests for four forward models, modified at 3xo from the best fit chemistry values inferred in the retrieval 
(on the low-resolution data). Top left: the metallicity is reduced to Z = 0.270. Top right: the metallicity is increased to Z = 
1.419. Bottom left: the C/O ratio is reduced to C/O = 0.658. Bottom right: the C/O ratio is increased to C/O = 0.898. In 
all four cases, the simulated phase-curve spectra do not match the observations, highlighting the sensitivity of those datasets to 
chemistry parameters. 
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D. COMPLEMENTARY FIGURES TO THE SECTION 5 


This appendix contains the complementary figures to the main article Section 5, Figure D1 to Figure D3. 
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Figure D1. ‘Transit (left) and eclipse (right) spectra of WASP-121 b analyzed in this work. Different observations are offset 
in the y-axis. Best-fit models from the 1D retrievals are shown in solid lines. Dashed lines show featureless models for visual 
comparison. 
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Figure D2. Posterior distributions obtained for the three transit fits. Yellow: Transit 1. Blue: Transit 2. Green: Transit 3. 
Note that the chemistry parameters of the equilibrium model GGCHEM in those fits are fixed to the median values obtained by 
the 1.5D phase-curve retrieval (log Z — -0.19, C/O — 0.80). Transit 1 shows an atmosphere with moderate hazes but absorption 
at high altitude from water. Transit 2 shows multimodal solutions involving either a fully ionized atmosphere with un-physically 
high temperatures, or a cloudy/hazy atmosphere. Transit 3 displays an atmosphere with opaque clouds. 
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Figure D3. Posterior distributions obtained for the five eclipse fits. Yellow: Eclipse 1. Blue: Eclipse 2. Green: Eclipse 3. Red: 
Eclipse 4. Purple: Eclipse 5. Note that the chemistry parameters of the equilibrium model GGCHEM in those fits are fixed to 
the median values obtained by the 1.5D phase-curve retrieval (log Z — -0.19, C/O — 0.80). The five eclipses are consistent with 
similarly inverted thermal structures, however the posterior distributions are not the same in the middle of the atmosphere (T1 
to T3). Large-scale atmospheric variability could creates those observable departures in the thermal profiles. 
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E. COMPLEMENTARY FIGURES TO THE SECTION 6 


This appendix contains the complementary figures for Section 6 of the main article, Figure E1 to Fig. EA. 
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Figure E1. Temperature maps of WASP-121b in Mollweide projection, obtained at p — 5x 10? Pa (top) and p — 10? Pa 
(bottom) for t € [40, 185] days. The figure is accompanied by two 1 min 23s videos, available online at the journal, showing 
the evolution of the atmosphere during 145 days in the simulations. The movie shows the highly time-variable atmosphere of 
WASP-121 b, expected from a high-resolution flow simulation. Note that the temperature ranges are different. 


Temperature (K) 


ATMOSPHERIC VARIABILITY OF WASP-121B 


Pressure = 1 bar 


T(K) [P= 1 bar, t =49 d] 


1625 1800 1975 


H2 [P=1 bar, t 249 d] 


H [P 1 bar, t 249 d] 


C — 


0.0000 — 0.0375 0.0750 0.1125 — 0.1500 


log(H~ ) [P= 1 bar, t =49 d] 


log(e-) [P= 1 bar, t 249 d] 


T (K) [P = 1 bar, t =62 d] 


1625 1800 1975 


H2 [P = 1 bar, t «62 d] 


H [P=1 bar, t 262 d] 


0.0375 0.0750 0.1125 0.1500 


log(H~ ) [P = 1 bar, t 262 d] 


log(e-) [P= 1 bar, t 262 d] 


Pressure = 50 mbar 


T (K) [P = 50 mbar, t =49 d] 


1700 2100 2500 


H2 [P 2 50 mbar, t =49 d] 


H [P = 50 mbar, t 249 d] 


0.0000 0.0375 0.0750 0.125 0.1500 


log(H ^ ) [P = 50 mbar, t 249 d] 


0.0000 


T (K) [P = 50 mbar, t =62 d] 


1700 2100 2500 


H2 [P 2 50 mbar, t 262 d] 


0.0375 0.0750 — 0.1125 


log(H~ ) [P = 50 mbar, t =62 d] 


-13.0 -12.0 -11.0 


log(e-) [P = 50 mbar, t =62 d] 


0.1500 


41 


Figure E2. WASP-121 b chemical maps centered around the sub-stellar point for Ho, H, H^ and e at two different times 
(t = 49 and t = 62 days) and two pressure levels (p = 10° Pa and p = 5x10? Pa). Those maps are obtained by post-processing 
the temperature fields (top row) from the BoB dynamics calculations with the TAUREX library. 
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Figure E3. WASP-121 b chemical maps centered around the sub-stellar point for H20, CO, CH4 and CO» at two different times 
(t = 49 and t = 62 days) and two pressure levels (p = 10° and p = 5 x 10? Pa). Those maps are obtained by post-processing 
the temperature fields (top row) obtained from the BoB dynamics calculations with the TAUREX library. 
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Figure E4. WASP-121 b chemical maps centered around the sub-stellar point for TiO, VO and FeH at two different times 
(t = 49 and t = 62 days) and two pressure levels (p = 10° and p = 5 x 10? Pa). Those maps are obtained by post-processing 
the temperature fields (top row) obtained from the BoB dynamics calculations with the TAUREX library. 


